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Abstract 

We simulate the natural frequencies and the acoustic wave propagation characteristics of graphene nanoribbons 
(GNRs) of the type (8,0) and (0,8) using an equivalent atomistic-continuum FE model previously developed by 
some of the authors, where the C-C bonds thickness and average equilibrium lengths during the dynamic loading 
are identified through the minimisation of the system Hamiltonian. A molecular mechanics model based on the 
UFF potential is used to benchmark the hybrid FE models developed. The acoustic wave dispersion characteristics 
of the GNRs are simulated using a Floquet-based wave technique used to predict the pass-stop bands of periodic 
mechanical structures. We show that the thickness and equilibrium lengths do depend on the specific vibration 
and dispersion mode considered, and that they are in general different from the classical constant values used in 
open literature (0.34 nm for thickness and 0.142 nm for equilibrium length). We also show the dependence of the 
wave dispersion characteristics versus the aspect ratio and edge configurations of the nanoribbons, with widening 
band-gaps that depend on the chirality of the configurations. The thickness, average equilibrium length and edge 
type have to be taken into account when nanoribbons are used to design nano-oscillators and novel types of 
mass sensors based on periodic arrangements of nanostructures. 
PACS 62.23.Kn • 62.25.Fg • 62.25.Jk 



Introduction 

Graphene nanoribbons (GNRs) [1] have attracted a sig- 
nificant interest in the nanoelectronics community as 
possible replacements to silicon semiconductors, quasi- 
THz oscillators and quantum dots [2], The electronic 
state of GNRs depend significantly on the edge struc- 
ture. The zigzag layout provides the edge localized state 
with non-bonding molecular orbitals near the Fermi 
energy, with induced large changes in optical and elec- 
tronic properties from quantization. DFT calculations 
and experimental measurements have shown that zigzag 
edge GNRs can show metallic or half-metallic behaviour 
(depending on the spin polarization in DFT simula- 
tions), while armchair nanoribbons are semiconducting 
with an energy gap decreasing with the increase of the 
GNR width [3-5]. GNRs have also been prototyped as 
photonics waveguides by Law et al. [6], and recently 
proposed for thermal phononics to control the reduc- 
tion of thermal conductivity by Yosevich and Savin [7], 
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In this study, we describe the mechanical vibration 
natural frequencies and acoustic wave dispersion char- 
acteristics of graphene nanoribbons considered as per- 
iodic structures. In structural dynamics design, the 
wave propagation characteristics of periodic systems 
(both ID and 2D) have been extensively used to tune 
the acoustic and vibrational signature of structures, 
materials and sensors [8-10], while at nanoscale level 
the periodicity of nanotubes array has also been used 
to develop nanophotonics crystals (see for example the 
study of Kempa and et al. [11]). Hod and Scuseria 
have also observed that the presence of a central 
mechanical load (or uniform inposed displacements) 
in bridged-bridged nanoribbons induces a significant 
electromechanical response in bending and torsional 
deformations [5]. We focus in this article on nanorib- 
bon architectures of the type (8,0) and (0,8). While the 
results present in this manuscript are related to these 
specific nanoribbon topologies, the general algorith 
that we proposed can be readily extended to analyse 
more general graphene architectures. The nanoribbon 
models are developed using a hybrid atomistic conti- 
nuum-Finite Element (FE) model (also called lattice 
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[12]), in which the carbon-carbon (C-C) covalent 
bonds are represented by Timoshenko structural 
beams with equivalent mechanical properties (Young's 
modulus and Poisson's ratio) derived by the minimisa- 
tion of the Hamiltonian of the structural system, or 
total potential energy for the static case [12-14]. It is 
worth to notice that the concept of the Hamiltonian of 
a system is not limited to problems associated to quan- 
tum mechanics, but it is also used in a large variety of 
variational problems related to the dynamics and stabi- 
lity of engineering and mechanical structures [15,16]. 
The equivalent mechanical properties for the sp 2 C-C 
bond are expressed in terms of the thickness of the 
bond itself. It is useful to reiterate that there is neither 
a physical thickness per se for the covalent bonds, nor 
for the carbon atoms involved in the bond. Nonethe- 
less, when subjected to a mechanical static loading, the 
nanostructure tends to reach its equilibrium state cor- 
responding to the minimum potential energy. The geo- 
metric and material configuration of the equivalent 
continuum mechanics structures used to represent the 
graphene (plates and/or shells) will be therefore be 
defined by the energy equilibrium conditions of the 
nanostructure, and cannot be ascribed as fixed. The 
length of the covalent bonds merits also some consid- 
erations. In finite size rectangular single layer gra- 
phene sheets (SLGS), the lengths of the C-C bonds at 
equilibrium after mechanical loading are unequal, 
ranging between 0.136 and 0.144 nm, and depend on 
the type of loading, size and boundary conditions 
[17,18], as well as the location on the SLGS itself (i.e. 
the edges [19]). This fact contrasts with the classical 
use of the fixed value of 0.142 nm at equilibrium con- 
sidered in most mechanical simulations [20-23]. The 
variation of the thickness and the distributions of 
lengths at equilibrium is important factors to consider 
when computing the homogenised mechanical proper- 
ties of the graphene, i.e. the equivalent mechanical 
performance of the graphene seen as a continuum. In 
this study, we will show that the thickness and the 
equilibrium length distributions assume some specific 
values in GNRs also when undergoing a mechanical 
resonant behaviour, both as a single nanostructure in 
free-free vibration conditions, and as periodic ele- 
ments in a one-dimensional (ID) acoustic wave pro- 
pagation case. However, the thickness and 
equilibrium lengths for the mechanical vibration case 
will be determined minimimsing the Hamiltonian of 
the system, rather that the total potential energy of 
the static loading case. Similar to the static in-plane 
and out-plane loading cases [12,13], those values can 
be different from the ones usually adopted in open 



literature. We will also show that the chirality of the 
GNRs (and their edge effects in nanoribbons with 
short widths) provides different acoustic wave disper- 
sion properties, which should be taken into account 
when GNRs are considered for potential nanoelectro- 
mechanical systems (NEMS) applications. 

Modeling 

Atomistic-FE model 

We use the atomistic-continuum equivalence model for 
the sp carbon-carbon bonds to extract the equivalent 
isotropic mechanical properties (Young's modulus and 
Poisson's ratio) as a functions of the thickness d of the 
C-C bond [13,14]. The model is based on the equiva- 
lence between the harmonic potential provided by force 
models such as AMBER or linearised Morse, and the 
strain energies associated to out-of-plane torsional, axial 
and bending deformation of a deep shear Timoshenko 
beam: 
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The first row of (1) corresponds to the equivalence 
between stretching and axial deformation mechanism 
(with E Y being the equivalent Young's modulus), while 
the second one equates the torsional deformation of the 
C-C bond with the pure shear deflection of the struc- 
tural beam associated to an equivalent shear modulus G. 
Contrary with analogous approaches previously used 
[21,23], the term related to the in-plane rotation of the 
C-C bond (third row of 1) is equated to a bending strain 
energy associated to a deep shear beam model, rather 
than a flexural one, to take into account the shear defor- 
mation of the cross section. The shear correction term 
becomes necessary when beams assume aspect ratios 
lower than 10 [24], which is the case for the C-C bonds 
with average lengths and thickness presented in in open 
literature (see the article of Huang et al. [25]). For circu- 
lar cross sections, the shear deformation constant can be 
expressed as [13]: 
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In (2), A s = A/F s is the reduced cross section of the 
beam by the shear correction term F s [26] : 



F s = 



6 + 12v + 6v 2 
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The insertion of (2) and (3) in (1) leads to an non- 
linear relation between the thickness d and the Poisson's 
ratio v of the equivalent beam [13]: 



k,d 2 4A + B 

kb = 

16 A + B 

where 

A = 112L 2 fe r + 192L 2 fc r v + 64L 2 fe T v 2 
B = 9k r d 2 + I8k r d 4 v + S)k r d 4 v 2 



(4) 
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The values for the force constants for the AMBER 
model are k r = 6.52 x 10" 7 N-mm" 1 , kg = 8.76 x 10' 10 N 
• nm • rad" 2 and k z = 2.78 x 10" 10 N • nm' 1 • rad' 2 . The 
equivalent mechanical properties of the C-C bond can 
be determined performing a nonlinear optimisation of 
(1) using a Marquardt algorithm. The C-C bond can 
then be discretised as a single two-nodes three-dimen- 
sional Finite Element model beam with a 6 x 6 stiffness 
matrix [K] e described in [27], where the nodes represent 
the atoms. The mass matrix [M] e of the bond is repre- 
sented through a lumped matrix approach [28]: 



[M] e = diag[^ 



m. 



0 0 0 



(7) 



3 3 3 

where m c = 1.9943 x 10' 26 kg. The elemental matrices 
are then assembled in the usual Finite Element fashion 
as global stiffness and mass matrices [K] and [M], 
respectively, which can be subsequently used to formu- 
late the undamped eigenvalue problem [29]: 



([K] 



? [M]){x} = {0} 



(8) 



Equation 8 is solved using a classical Block Lanczos algo- 
rithm implemented in the commercial FE code ANSYS 
(Rel. 12). According to Equation 2-4, the natural frequen- 
cies £0, are, however, dependent on the thickness d. In the 
hybrid FE simulation, we consider also the variation of the 
average bond length / across the graphene sheet, a phenom- 
enon observed in several models of SLGSs subjected to 
mechanical loading [13,17,19,30]. To identify a unique set 
of thickness and equilibrium lengths for a specific eigenso- 
lution, we minimise the Hamiltonian of the system [15]: 



H = T+U 



(9) 



where T and U are the kinetic and strain energies of 
the system, respectively. Using the mass-normalized 
normal modes [O] associated to the eigenvalue problem 
[29], the Hamiltonian (9) for each eigensolution i can be 
rewritten as: 



The ID wave propagation analysis is carried out using 
a technique implemented by Tee et al. [10] and Aberg 
and Gudmundson [31]. Applying the Floquet conditions 
between the left and the right nodal degrees of freedom 
(DOFs) {u} L and {uf one obtains: 



e- ik '{uf 



(11) 



where -n < k x < tc is the propagation constant within 
the first Brillouin zone [32]. The generalized DOFs of 
the system will be complex (real and imaginary part), 
while for traveling waves the propagation constant k x 
will be solely real [32]. Equation 11 can be, therefore, 
recast as: 



Ullm = {Ulfm cosk x - { u }Re sin k x 
(Re = {Ulfm cos k * + MlL sin k x 



(12) 



The real and imaginary parts of the domain in the FE 
representation are produced creating two superimposed 
meshes, linked by the boundary conditions [10,31] (12). 
For a given wave propagation constant k x , the resultant 
eigenvalue problem provides the frequency associated to 
the acoustic wave dispersion curve. Similar to the 
undamped eigenvalue problem, the minimisation of the 
Hamiltonian (10) is also carried out for the wave propa- 
gation case to identify the set of thickness and average 
bond length required for the eigenvalue solution. 

Molecular mechanics approach 

The molecular mechanics (MM) simulations were per- 
formed with Gaussian [33], using the universal force 
field (UFF) developed by Rappe et al. [34]. Force-field- 
based simulations are convenient to represent the acous- 
tic/mechanical dynamics behaviour, because they use 
explicit expressions for the potential energy surface of a 
molecule as a function of the atomic coordinates. The 
UFF is also well suited for dynamics simulations, allow- 
ing more accurate vibration measurements than many 
other force fields, which do not distinguish bond 
strengths. The UFF is a purely harmonic force field with 
a potential-energy expression of the form: 

E = ^2 e r + ^2 e 6 +^2 E <t> +^2 E u £vdw + XI £ei ( 13 ' 

The valence interactions consist of bond stretching 
(E R ), which is a harmonic term and angular distortions. 
The angular distortions are the bond angle bending (Eg), 
described by a three-term Fourier cosine expansion, the 
dihedral angle torsion (E^) and inversion terms (out-of- 
plane bending) (E cu ). Ey and E cu are described by cosine- 
Fourier expansion terms. The non-bonded interactions 
consist of van der Waals (£vdw) and electrostatic (£ e i) 
terms. £vdw are described by a Lennard-Jones potential, 
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while E el described by a Coulombic term. The functional 
form of the above energy terms is given as follows: 



Er 
Eg 

C 2 

Ci 
Co 

£vdw 



fei(r-r 0 f 

k 2 (Co + C\ cos 9 + d cos 20) 
1 



4 sin 2 0 



-4C 2 cos0 
C 2 (2cos 2 0 + 1) 

ks(l ± COSH0) 

1 ±cos(n0)] 

12 

- 2 



(14) 



fe 4 
D 



gig/ 



Here ki, ki, k^ and /c 4 are force constants, 6>o is the 
natural bond angle, D is the van der Waals well depth, 
r" is the van der Waals length, q t is the net charge of an 
atom, e is the dielectric constant and r» is the distance 
between two atoms. In nanotubes, the atoms have no 
net charge, so the £ e i term is always zero. The torsion 
term, E v , turns out to be of great importance. Detailed 
values of these parameters in Equation 14 can be found 
in Ref. [34]. Some of the authors have successfully used 



a similar MM approach to describe the mechanical 
vibrations of single-walled carbon nanotubes [35] and 
boron-nitride nanotubes [36]. Other molecular 
mechanics approaches have been successfully used to 
describe the structural mechanics aspects of SWCNTs 
and MWCNTs (see for example Sears and Batra [37]). 

Results and discussions 

Molecular mechanics and atomistic-FE models 

Figure 1 shows the comparison between the MM simu- 
lations and the results from the hybrid FE models for a 
(8,0) nanoribbon at different lengths (6.03, 12.18, 18.34 
and 24.49 nm). The equilibrium lengths are / = 0.142 
nm for all cases considered. For the flexural modes the 
hybrid FE approach identifies a bond thickness d of 
0.077 nm, with only a 3% difference from the analogous 
thickness value assocoated to the first torsional mode is 
considered. The identified thickness value compares well 
with the 0.074-0.099 nm found by some of the authors 
in uni-axial tensile loading cases related to single layer 
graphene sheets [13], with the 0.0734 nm in uni-axial 
stretching using first generation Brenner potential [25], 
and the 0.0894 nm identified by Kudin et al using ab 
initio techniques [38]. Gupta and Batra [39] find a 
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Figure 1 Comparison between MM (full markers) and hybrid-FE (empty markers) natural frequencies for (8,0) SLGSs with different 
widths. 
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thickness of 0.080 nm for the con frequency of a fully 
clamped single layer graphene sheet (SLGS) with dimen- 
sions 3.23 nm x 2.18 nm, combining a MD simulation 
and results from the continuum elasticity of plates. It is 
worth to notice that these results are significantly differ- 
ent from the usual 0.34 nm inter-atomic layer distance 
adopted by the vast majority of the research community 
in nanomechanical simulations. The percentage differ- 
ence between our MM and hybrid FE natural frequen- 
cies is on average around 3 for all the flexural modes. 
The torsional frequencies for the nanoribbons with the 
lowest aspect ratio provide a higher error (5%), suggest- 
ing that the assumption of equal in-plane and out-of- 
plane torsional stiffness with the AMBER model in 
Equation 1 leads to a slightly lower out-of-plane tor- 
sional stiffness of the nanoribbon. 

Wave propagation in bridged nanoribbons with different 
chirality 

The ID wave propagation analysis has been carried out 
on (8,0) nanoribbons with a length of 15.854 nm along 
the zigzag direction, and 15.407 nm along the armchair 
direction for the (0,8) cases. The hybrid FE models have 
been subjected to simply supported (SS) conditions, 
clamping the relevant DOFs in the middle location of 
the ribbons, and allowing, therefore, to apply the rela- 
tions (12) using a set of constraint equations. The wave 
dispersion characteristics for the propagation along the 
zigzag edge of the nanoribbons for the first Brillouin 
zone [32] are shown in Figure 2. The mode shapes asso- 
ciated to the first four pass-stop bands (Figure 3) are 
typical of periodic SS structural beams under bending 
deformation [40], while from our observations the out- 
of-plane torsional modes appear for the 5th and 6th 
wave dispersion characteristics. 

A more significant discrepancy between wave disper- 
sion curves can be observed in Figure 2, when compar- 
ing the pass-stop band behaviour for the propagation 
along the zigzag and armchair directions. Only the first 
acoustic flexural wave dispersion characteristic is vir- 
tually unchanged, while for the other curves we observe 
a strong decrease in terms of magnitude, as well as 
mode inversion. The first stop band is significantly 
decreased by 25 GHz for k x = n - the armchair case 
gives a frequency drop of 39 GHz for the same propaga- 
tion constant. Similar decreases in band gaps are 
observed for higher frequencies, while mode inversion 
(flexural to torsional) is observed for the armchair pro- 
pagation around kjn = 0.42, while for the armchair 
case the mode inversion is located around 0.8 kjn. 
From the mechanical point of view, a possible explana- 
tion for this peculiar behaviour can be given considering 
the intrinsic anisotropy of the in-plane properties of 
finite size graphene sheets. Reddy et al. [17,41] have 



observed anisotropy ratios between 0.92 and 0.94 in 
almost square graphene sheets subjected to uni-axial 
loading, while similar orthotropic ratios have been iden- 
tified also by Scarpa et al. [13]. The GNRs considered 
here have an aspect ratio close to 6, which induces the 
edges to provide a higher contribution to the homoge- 
nized mechanical properties due to Saint Venant effects 
[42]. A further confirmation of the effective in-plane 
mechanical anisotropy on the GNRs is apparent also 
from the non-dimensional dispersion curves shown in 
Figure 4. For that specific case, the GNRs have one side 
fixed (1.598 nm for the armchair, and 1.349 nm for the 
zigzag), with minimized thickness d equal to 0.074 and 
0.077 nm and C-C bond equilibrium lengths of / = 
0.142 nm for the armchair and zigzag cases, respectively. 
The dimensions of the nanoribbons are varied adjusting 
the aspect ratios (2.4 and 8), to obtain armchair and zig- 
zag GNRs with similar dimensions. We have further 
nondimensionalised the dispersion curves using the 
values of the first dispersion relation (m 0 ) for the arm- 
chair configuration at k x = n/4. The GNR with an aspect 
ratio of 2.4 (Figure 4a) shows significant difference s in 
terms of dispersion characteristics between the armchair 
and the zigzag configurations, with a reduced band-gap 
of A(co/oj 0 ) equal to 3 for the armchair, against the value 
of 5 for the zigzag at the end of the first Brillouin zone 
{kjn = 1). Between 4 <co/co 0 <10, the wave dispersions 
appear to be composed by combinations of flexural 
plate-like modes with torsional components, with mode 
veering occurring between 0.45 < kjn < 0.65. The zig- 
zag-edged GNRs tend to show a narrowing of the non- 
dimensional dispersion characteristics within the same 
(o/a>o range considered. At higher non-dimensional wave 
dispersions, both armchair and zigzag nanoribbons tend 
to show beam-like dispersion characteristics [8,40,43]. 
The nanoribbons with higher aspect ratio (Figure 4b) 
show the pass-stop band behaviour typical of SS peri- 
odic structures made of Euler-Bernoulli beams [40]. 
However, while the first non-dimensional dispersion 
curve is identical, the following dispersion characteristics 
show a marked difference between zigzag and armchair 
configurations, with the zigzag GNRs having the highest 
co/co 0 values. It is also worth of notice that while the zig- 
zag configuration shows a dispersion curve provided by 
a torsional wave (straight line between 0 < kjn < 0.62 
at m/coo = 37.4), the armchair GNR appears to be gov- 
erned by flexural waves within the non-dimensional fre- 
quency interval considered. This type of behaviour 
suggests also that the specific morphology of the edges 
(combined with the small transversal dimensions of the 
GNRs) affect the acoustic wave propagation characteris- 
tics, both contributing to an overall mechanical aniso- 
tropy of the equivalent beams, as well as providing 
specific wave dispersion characteristics at higher 
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(a) 



(b) 





(c) 



(d) 



Figure 3 Mode shapes (real parts) for a (8,0) GNR (length 15.854 nm) with propagation constant k x = n/4 along the zigzag direction. 

(a) co, = 8.84 GHz; (b) co 2 = 29.35 GHz; (c) fo 3 = 52.5 GHz; (d) cu 4 = 82.5 GHz. 



frequencies. Moreover, the widening of the band gap 
observed in Figure 2 for the armchair configuration 
recalls some similarity to the variation of the energy gap 
of the electronic states noted in analogous armchair 
GNRs [3]. For a fixed width of 2.25 nm and and aspect 
ratio of 2.4 (i.e. 5.4 nm), the first pass-stop band at k x = 
0 is located at 180 GHz. For the same fixed width but 
higher aspect ratio (8.0, corresponding to a transverse 
length of 18 nm), the same pass-stop band first fre- 
quency for k x = 0 is equal to 15 GHz, 12 times lower 
than the low aspect ratio case (Figure 4). Moreover, for 
the higher aspect ratio we observe aAcv = 18 GHz, while 
the lower aspect ratio provides a pass-stop band fre- 
quency interval Aco = 90 GHz, five times higher when 
compared for the armchair nanoribbons at AR = 2.4. 
Passing between lengths of 0.25 and 3 nm, Barone, Hod 
and Scuseria observe a decrease in energy gab by a fac- 
tor of 3 for bare PBEs, and by 5 for bare HSEs [3]. 
When we consider the variation of the energy of the 



system proportional to the kinetic energy (and therefore 
approximately Am 2 ), ther ratio of the pass-stop bands 
for the armchair nanoribbons with different aspect ratios 
is compatible with the decrease of energy gap observed 
through DFT simulations [3]. 

Conclusions 

In this study, we have presented a new methodology to 
derive the mechanical structural dynamics characteris- 
tics and acoustic wave dispersion relations for graphene 
nanoribbons using an hybrid Finite Element approach. 
The technique, benchmarked against a Molecular 
Mechanics model, allows to identify the mechanical nat- 
ural frequencies and associated modes shapes, as well as 
the pass-stop band acoustic characteristics of periodic 
arrays of GNRs. 

The numerical results from the minimisation of the 
Hamiltonian in the hybrid FE method show that the 
commonly used value in nanomechanical simulations 
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Figure 4 Non-dimensional dispersion curves for zigzag (continuous lines) and armchair (dashed lines) (8,0) GNRs with different aspect 
ratios (AR). All the results are minimized for the thickness d only. 



for the thickness (0.34 nm) is not adequate to represent 
the effective structural dynamics of the system. Thick- 
ness values identified through the minimisation of the 
Hamiltonian vary in a restricted range around 0.07 nm 
for the AMBER force model used in this study. We also 



observe a distribution of the C-C bond lengths corre- 
sponding to average values between 0.142 nm and 0.145 
nm, after the minimisation for specific modes. However, 
the minimised thickness does not show any particular 
dependence over the type of mode shape considered. 
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Only for pure torsional modes a small percentage varia- 
tion from the baseline d = 0.074 nm value is observed. 

We also show that graphene manoribbons exhibit a 
significant dependence of the acoustic wave propagation 
properties over the type of edge and aspect ratio, quite 
similarly to what observed for their electronic state. 
This feature suggests a possible combined electro- 
mechanical approach to design multifunctional wave- 
guide-type band filters. 

The use of periodic assemblies of graphene nanorib- 
bons seems also a design feature that could lead to 
potential breakthroughs in terms of mass-sensors con- 
cepts, with enhanced selectivity provided by the periodic 
distribution of constraints and supports. The model pro- 
posed in this study allows to design and simulate these 
novel devices. 
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